library(dplyr)
library(lavaan)
library(DiagrammeR)
library(tidyr)
library(ggplot2)

# For saving SEM diagrams:
library(purrr)
library(DiagrammeRsvg)
library(rsvg)
library(png)
library(grid)
library(ggpubr)

Import data

combined=read.csv("data/annual_averages/annual_data_compiled_regions.csv",stringsAsFactors = F)
cnames=read.csv("analysis/column_names_region.csv", stringsAsFactors = F)
dsub=filter(combined, Year>=1975) %>% arrange(Region,Year)
focaldata=dsub[,cnames$Datacolumn]
fvars=cnames$Shortname
colnames(focaldata)=fvars
regions=unique(focaldata$region)
regionorder=c("West","North","South")
years=1975:2021

focaldata = focaldata %>% 
  mutate(tzoop=hcope+clad+mysid+pcope+rotif_m,
         tzoop_e=hcope_e+clad_e+mysid_e+pcope_e+rotif_e,
         hzoop=hcope+clad+rotif_m,
         hzoop_e=hcope_e+clad_e+rotif_e,
         pzoop=mysid+pcope,
         pzoop_e=mysid_e+pcope_e,
         turbid=-secchi) 
fvars=c(fvars,"tzoop","tzoop_e",
        "hzoop","hzoop_e",
        "pzoop","pzoop_e","turbid")
cnames=rbind(cnames,data.frame(Longname=NA,Shortname=c("tzoop","tzoop_e",
                                                       "hzoop","hzoop_e",
                                                       "pzoop","pzoop_e","turbid"),
                               Diagramname=c("total zooplankton",
                                             "total zooplankton\nenergy",
                                             "herbivorous\nzooplankton",
                                             "herbivorous\nzooplankton\nenergy",
                                             "predatory\nzooplankton",
                                             "predatory\nzooplankton\nenergy",
                                             "turbidity"),
                               Datacolumn=NA,Log=c(rep("yes",6),"no"),
                               Color=c("black","black","#ED7D31","#ED7D31","#7030A0",
                                       "#7030A0","#4472C4")))

#focal variables
varnames=c("temp","flow","turbid","chla","hzoop","pzoop","tzoop","amphi","potam","corbic","estfish","estfish_bsmt","estfish_stn")

source("analysis/myLavaanPlot.r")
source("analysis/semDiagramFunctions.r")

Data prep

Log transform, scale

#log transform
logvars=fvars[cnames$Log=="yes"]
logtrans=function(x) {
  x2=x[which(!is.na(x))]
  if(any(x2==0)) {log(x+min(x2[which(x2>0)],na.rm=T))}
  else {log(x)}
}
focaldatalog = focaldata %>% 
  mutate(flow=flow-min(flow,na.rm=T)) %>%  #get rid of negative flow values
  mutate_at(logvars,logtrans)

#scale data
fdr0=focaldatalog
tvars=fvars[-(1:2)]

fdr=fdr0 %>% group_by(region) %>% 
  #lag
  #mutate_at(tvars,list("1"=lag)) %>% 
  #scale
  mutate_at(-(1:2),scale) %>% 
  ungroup() %>% 
  as.data.frame() 

#detrended data
fdr_dtr=fdr0 %>% group_by(region) %>% 
  #detrend
  mutate_at(tvars,function(x) { 
    x<<-x
    if(!all(is.na(x))) {
      if((length(which(x==0))/length(x))<0.5) {
        x2<<-x
        x2[x2==0]=NA
        res<<-residuals(lm(x2~years))
        out=x
        out[which(!is.na(x2))]=res
        return(out)
      } else {return(x)}
    } else {return(x)}
  }) %>%
  #lag
  #mutate_at(tvars,list("1"=lag)) %>% 
  #scale
  mutate_at(-(1:2),scale) %>% 
  ungroup() %>% 
  as.data.frame() 

Time series plots

## Note: Using an external vector in selections is ambiguous.
## i Use `all_of(varnames)` instead of `varnames` to silence this message.
## i See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.

Other useful plots

Breakdown of total zooplankton biomass.

## Warning: Removed 15 rows containing missing values (position_stack).

Similarity of fish indices.

## Warning: Removed 7 row(s) containing missing values (geom_path).

## Warning: Removed 7 row(s) containing missing values (geom_path).

Correlation between biomass and energy.

for(i in 1:length(regions)) {
  dtemp=filter(fdr,region==regions[i])
  print(regions[i])
  print(cor(dtemp$tzoop,dtemp$tzoop_e,use = "p"))
  print(cor(dtemp$hzoop,dtemp$hzoop_e,use = "p"))
  print(cor(dtemp$pzoop,dtemp$pzoop_e,use = "p"))
}
## [1] "North"
##           [,1]
## [1,] 0.9989681
##           [,1]
## [1,] 0.9958644
##           [,1]
## [1,] 0.9990414
## [1] "South"
##           [,1]
## [1,] 0.9987488
##           [,1]
## [1,] 0.9977845
##           [,1]
## [1,] 0.9993711
## [1] "West"
##          [,1]
## [1,] 0.999451
##           [,1]
## [1,] 0.9986605
##           [,1]
## [1,] 0.9994171

Cross-correlation matrices

(only sig correlations shown… no correction for multiple comparisons)

## Warning: `expand_scale()` is deprecated; use `expansion()` instead.

## Warning: `expand_scale()` is deprecated; use `expansion()` instead.

## Warning: `expand_scale()` is deprecated; use `expansion()` instead.

## Warning: `expand_scale()` is deprecated; use `expansion()` instead.

Note correlation of fish indices, or lack thereof.

SEM model

With and without detrending.

West

#1
# modwest='zoop=~hcope+mysid
#         fish=~estfish_bsmt+estfish_bsot
#         zoop~chla+potam+flow
#         chla~potam+flow
#         fish~zoop+flow
# '
#2
# modwest='chla~potam+flow
#         tzoop~chla+potam+flow
#         estfish_bsmt~tzoop+flow
#         estfish_bsot~tzoop+flow
# '
#3
# modwest='chla~potam+flow+temp+turbid
#         tzoop~chla+potam+flow+temp+turbid
#         amphi~chla+potam+flow+temp+turbid
#         estfish_bsmt~tzoop+amphi+flow
#         #estfish_bsot~tzoop+amphi+flow+temp+turbid
#         amphi~~tzoop
# '
#4
# modwest='chla~potam+flow+temp+turbid
#         hzoop~chla+potam+flow+temp+turbid
#         pzoop~chla+potam+flow+temp+turbid+hzoop
#         amphi~chla+potam+flow+temp+turbid
#         estfish_bsmt~hzoop+pzoop+amphi+flow+temp+turbid
#         amphi~~hzoop+pzoop
# '
#5
modwest='chla~potam+flow+temp+turbid
        hzoop~chla+potam+flow+temp+turbid
        pzoop~chla+potam+flow+temp+turbid+hzoop
        fish~hzoop+pzoop+potam+flow+temp+turbid
        fish=~estfish+estfish_stn+estfish_bsmt
'

modfitwest=sem(modwest, data=filter(fdr,region=="West"))
modfitwest_dtr=sem(modwest, data=filter(fdr_dtr,region=="West"))
summary(modfitwest, standardized=T, rsq=T)
## lavaan 0.6-10 ended normally after 40 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        30
##                                                       
##                                                   Used       Total
##   Number of observations                            40          47
##                                                                   
## Model Test User Model:
##                                                       
##   Test statistic                                21.159
##   Degrees of freedom                                15
##   P-value (Chi-square)                           0.132
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fish =~                                                               
##     estfish           1.000                               0.875    0.879
##     estfish_stn       0.919    0.115    8.016    0.000    0.804    0.883
##     estfish_bsmt      0.848    0.143    5.943    0.000    0.742    0.752
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   chla ~                                                                
##     potam            -0.285    0.169   -1.688    0.091   -0.285   -0.285
##     flow              0.281    0.182    1.540    0.124    0.281    0.280
##     temp             -0.281    0.150   -1.870    0.062   -0.281   -0.278
##     turbid           -0.221    0.199   -1.112    0.266   -0.221   -0.213
##   hzoop ~                                                               
##     chla              0.520    0.103    5.067    0.000    0.520    0.599
##     potam            -0.277    0.114   -2.439    0.015   -0.277   -0.319
##     flow              0.137    0.122    1.120    0.263    0.137    0.157
##     temp              0.025    0.102    0.245    0.807    0.025    0.028
##     turbid           -0.220    0.131   -1.682    0.093   -0.220   -0.245
##   pzoop ~                                                               
##     chla              0.397    0.084    4.727    0.000    0.397    0.455
##     potam            -0.034    0.078   -0.437    0.662   -0.034   -0.039
##     flow             -0.344    0.079   -4.351    0.000   -0.344   -0.392
##     temp             -0.005    0.065   -0.075    0.940   -0.005   -0.006
##     turbid            0.109    0.087    1.262    0.207    0.109    0.121
##     hzoop             0.594    0.101    5.889    0.000    0.594    0.591
##   fish ~                                                                
##     hzoop             0.214    0.157    1.364    0.173    0.244    0.212
##     pzoop             0.299    0.145    2.068    0.039    0.342    0.298
##     potam            -0.301    0.090   -3.332    0.001   -0.343   -0.343
##     flow              0.195    0.101    1.928    0.054    0.222    0.221
##     temp              0.010    0.072    0.139    0.890    0.011    0.011
##     turbid            0.221    0.101    2.202    0.028    0.253    0.244
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .estfish           0.225    0.064    3.517    0.000    0.225    0.227
##    .estfish_stn       0.182    0.053    3.469    0.001    0.182    0.220
##    .estfish_bsmt      0.424    0.102    4.147    0.000    0.424    0.435
##    .chla              0.718    0.161    4.472    0.000    0.718    0.719
##    .hzoop             0.303    0.068    4.472    0.000    0.303    0.402
##    .pzoop             0.123    0.028    4.472    0.000    0.123    0.162
##    .fish              0.065    0.040    1.624    0.104    0.085    0.085
## 
## R-Square:
##                    Estimate
##     estfish           0.773
##     estfish_stn       0.780
##     estfish_bsmt      0.565
##     chla              0.281
##     hzoop             0.598
##     pzoop             0.838
##     fish              0.915
summary(modfitwest_dtr, standardized=T, rsq=T)
## lavaan 0.6-10 ended normally after 33 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        30
##                                                       
##                                                   Used       Total
##   Number of observations                            40          47
##                                                                   
## Model Test User Model:
##                                                       
##   Test statistic                                15.035
##   Degrees of freedom                                15
##   P-value (Chi-square)                           0.449
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fish =~                                                               
##     estfish           1.000                               0.648    0.702
##     estfish_stn       1.193    0.249    4.787    0.000    0.773    0.784
##     estfish_bsmt      0.989    0.249    3.972    0.000    0.640    0.648
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   chla ~                                                                
##     potam             0.175    0.157    1.119    0.263    0.175    0.178
##     flow              0.632    0.196    3.228    0.001    0.632    0.629
##     temp             -0.330    0.144   -2.295    0.022   -0.330   -0.325
##     turbid           -0.422    0.183   -2.307    0.021   -0.422   -0.418
##   hzoop ~                                                               
##     chla              0.465    0.133    3.497    0.000    0.465    0.465
##     potam            -0.013    0.134   -0.099    0.921   -0.013   -0.013
##     flow              0.466    0.185    2.521    0.012    0.466    0.464
##     temp             -0.026    0.128   -0.200    0.842   -0.026   -0.025
##     turbid           -0.460    0.164   -2.809    0.005   -0.460   -0.455
##   pzoop ~                                                               
##     chla              0.431    0.095    4.531    0.000    0.431    0.511
##     potam             0.019    0.084    0.230    0.818    0.019    0.023
##     flow             -0.269    0.124   -2.161    0.031   -0.269   -0.318
##     temp             -0.033    0.080   -0.406    0.685   -0.033   -0.038
##     turbid           -0.004    0.112   -0.033    0.973   -0.004   -0.004
##     hzoop             0.418    0.099    4.223    0.000    0.418    0.495
##   fish ~                                                                
##     hzoop             0.039    0.105    0.371    0.710    0.060    0.060
##     pzoop             0.227    0.118    1.934    0.053    0.351    0.295
##     potam            -0.085    0.074   -1.147    0.251   -0.132   -0.133
##     flow              0.524    0.131    3.998    0.000    0.809    0.803
##     temp             -0.041    0.070   -0.587    0.557   -0.063   -0.062
##     turbid           -0.018    0.099   -0.179    0.858   -0.027   -0.027
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .estfish           0.431    0.106    4.073    0.000    0.431    0.507
##    .estfish_stn       0.373    0.102    3.652    0.000    0.373    0.385
##    .estfish_bsmt      0.565    0.134    4.209    0.000    0.565    0.580
##    .chla              0.698    0.156    4.472    0.000    0.698    0.703
##    .hzoop             0.494    0.110    4.472    0.000    0.494    0.497
##    .pzoop             0.193    0.043    4.472    0.000    0.193    0.274
##    .fish              0.023    0.041    0.563    0.573    0.055    0.055
## 
## R-Square:
##                    Estimate
##     estfish           0.493
##     estfish_stn       0.615
##     estfish_bsmt      0.420
##     chla              0.297
##     hzoop             0.503
##     pzoop             0.726
##     fish              0.945
labelswest <- createLabels(modfitwest, cnames)

# residuals(modfitwest)
# modificationindices(modfitwest)

North

#1
# modnorth='zoop=~hcope+mysid
#         #fish=~estfish_bsmt+estfish_bsot
#         zoop~chla+potam+flow
#         chla~potam+flow
#         estfish_bsmt~zoop+flow
#         estfish_bsot~zoop+flow
# '
# modnorth='zoop=~clad
#         zoop~chla+corbic+potam+flow
#         chla~corbic+potam+flow
#         estfish_bsmt~zoop+flow+sside+chla
#         estfish_bsot~zoop+flow+sside+chla
# '
#2
# modnorth='chla~corbic+potam+flow
#         tzoop~chla+corbic+potam+flow
#         estfish_bsmt~tzoop+flow+chla
#         estfish_bsot~tzoop+flow
# '
#3
# modnorth='chla~corbic+potam+flow+temp+turbid
#         tzoop~chla+corbic+potam+flow+temp+turbid
#         amphi~chla+corbic+potam+flow+temp+turbid
#         estfish_bsmt~tzoop+amphi+flow+temp+turbid+chla+sside
#         #estfish_bsot~tzoop+amphi+flow+temp+turbid+sside
#         amphi~~tzoop
# '
#4
# modnorth='chla~corbic+flow+temp+turbid
#         hzoop~chla+corbic+flow+temp+turbid
#         pzoop~chla+corbic+flow+temp+turbid+hzoop
#         amphi~chla+corbic+flow+temp+turbid
#         estfish_bsmt~hzoop+pzoop+amphi+flow+temp+turbid+chla+sside
#         amphi~~hzoop+pzoop
# '
#5
modnorth='chla~corbic+flow+temp+turbid
        hzoop~chla+corbic+flow+temp+turbid
        pzoop~chla+corbic+flow+temp+turbid+hzoop
        fish~chla+hzoop+pzoop+corbic+flow+temp+turbid
        fish=~estfish+estfish_stn #+estfish_bsmt
        estfish_stn~~chla
'

modfitnorth=sem(modnorth, data=filter(fdr,region=="North"))
modfitnorth_dtr=sem(modnorth, data=filter(fdr_dtr,region=="North"))
summary(modfitnorth, standardized=T, rsq=T)
## lavaan 0.6-10 ended normally after 34 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        30
##                                                       
##                                                   Used       Total
##   Number of observations                            41          47
##                                                                   
## Model Test User Model:
##                                                       
##   Test statistic                                 4.919
##   Degrees of freedom                                 5
##   P-value (Chi-square)                           0.426
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fish =~                                                               
##     estfish           1.000                               0.833    0.812
##     estfish_stn       0.892    0.171    5.207    0.000    0.743    0.839
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   chla ~                                                                
##     corbic            0.449    0.126    3.549    0.000    0.449    0.450
##     flow              0.211    0.146    1.441    0.150    0.211    0.201
##     temp              0.182    0.131    1.392    0.164    0.182    0.178
##     turbid            0.106    0.134    0.793    0.428    0.106    0.103
##   hzoop ~                                                               
##     chla              0.131    0.120    1.095    0.273    0.131    0.152
##     corbic            0.316    0.120    2.628    0.009    0.316    0.367
##     flow             -0.054    0.127   -0.424    0.672   -0.054   -0.060
##     temp              0.033    0.113    0.288    0.773    0.033    0.037
##     turbid            0.402    0.108    3.713    0.000    0.402    0.454
##   pzoop ~                                                               
##     chla              0.580    0.089    6.508    0.000    0.580    0.602
##     corbic            0.346    0.095    3.635    0.000    0.346    0.360
##     flow             -0.496    0.093   -5.332    0.000   -0.496   -0.492
##     temp              0.003    0.083    0.031    0.975    0.003    0.003
##     turbid            0.122    0.092    1.335    0.182    0.122    0.124
##     hzoop             0.160    0.114    1.397    0.162    0.160    0.143
##   fish ~                                                                
##     chla             -0.554    0.187   -2.965    0.003   -0.666   -0.675
##     hzoop             0.460    0.157    2.925    0.003    0.553    0.483
##     pzoop             0.526    0.207    2.541    0.011    0.631    0.617
##     corbic           -0.266    0.153   -1.743    0.081   -0.320   -0.325
##     flow              0.232    0.161    1.443    0.149    0.279    0.271
##     temp              0.104    0.112    0.924    0.356    0.124    0.123
##     turbid            0.285    0.125    2.281    0.023    0.342    0.338
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##  .estfish_stn ~~                                                        
##    .chla              0.280    0.113    2.477    0.013    0.280    0.589
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .estfish           0.357    0.115    3.096    0.002    0.357    0.340
##    .estfish_stn       0.321    0.109    2.955    0.003    0.321    0.410
##    .chla              0.706    0.156    4.528    0.000    0.706    0.686
##    .hzoop             0.417    0.092    4.528    0.000    0.417    0.545
##    .pzoop             0.224    0.049    4.528    0.000    0.224    0.234
##    .fish              0.210    0.096    2.181    0.029    0.302    0.302
## 
## R-Square:
##                    Estimate
##     estfish           0.660
##     estfish_stn       0.590
##     chla              0.314
##     hzoop             0.455
##     pzoop             0.766
##     fish              0.698
summary(modfitnorth_dtr, standardized=T, rsq=T)
## lavaan 0.6-10 ended normally after 34 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        30
##                                                       
##                                                   Used       Total
##   Number of observations                            41          47
##                                                                   
## Model Test User Model:
##                                                       
##   Test statistic                                 5.852
##   Degrees of freedom                                 5
##   P-value (Chi-square)                           0.321
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fish =~                                                               
##     estfish           1.000                               0.608    0.613
##     estfish_stn       1.198    0.430    2.784    0.005    0.729    0.729
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   chla ~                                                                
##     corbic            0.515    0.158    3.258    0.001    0.515    0.485
##     flow              0.151    0.160    0.946    0.344    0.151    0.144
##     temp              0.173    0.134    1.290    0.197    0.173    0.167
##     turbid            0.200    0.142    1.412    0.158    0.200    0.193
##   hzoop ~                                                               
##     chla              0.342    0.166    2.063    0.039    0.342    0.350
##     corbic           -0.027    0.196   -0.138    0.890   -0.027   -0.026
##     flow              0.151    0.190    0.795    0.426    0.151    0.146
##     temp              0.091    0.159    0.569    0.569    0.091    0.089
##     turbid           -0.008    0.170   -0.046    0.963   -0.008   -0.008
##   pzoop ~                                                               
##     chla              0.745    0.107    6.955    0.000    0.745    0.725
##     corbic            0.292    0.121    2.418    0.016    0.292    0.267
##     flow             -0.544    0.117   -4.637    0.000   -0.544   -0.501
##     temp              0.019    0.098    0.198    0.843    0.019    0.018
##     turbid            0.006    0.105    0.056    0.955    0.006    0.006
##     hzoop             0.037    0.096    0.382    0.703    0.037    0.035
##   fish ~                                                                
##     chla             -0.411    0.215   -1.911    0.056   -0.676   -0.692
##     hzoop             0.073    0.101    0.724    0.469    0.120    0.120
##     pzoop             0.328    0.182    1.800    0.072    0.539    0.568
##     corbic           -0.412    0.190   -2.172    0.030   -0.678   -0.653
##     flow              0.326    0.177    1.840    0.066    0.537    0.521
##     temp              0.144    0.118    1.227    0.220    0.237    0.234
##     turbid            0.046    0.123    0.377    0.706    0.076    0.075
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##  .estfish_stn ~~                                                        
##    .chla              0.461    0.185    2.486    0.013    0.461    0.679
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .estfish           0.614    0.170    3.609    0.000    0.614    0.624
##    .estfish_stn       0.620    0.242    2.560    0.010    0.620    0.621
##    .chla              0.743    0.164    4.528    0.000    0.743    0.708
##    .hzoop             0.839    0.185    4.528    0.000    0.839    0.835
##    .pzoop             0.317    0.070    4.528    0.000    0.317    0.286
##    .fish              0.167    0.109    1.533    0.125    0.453    0.453
## 
## R-Square:
##                    Estimate
##     estfish           0.376
##     estfish_stn       0.379
##     chla              0.292
##     hzoop             0.165
##     pzoop             0.714
##     fish              0.547
labelsnorth <- createLabels(modfitnorth, cnames)

# residuals(modfitnorth)
# modificationindices(modfitnorth)

South

#1
# modsouth='zoop=~hcope+mysid
#         #fish=~estfish_bsmt+estfish_bsot
#         zoop~chla+corbic+flow
#         chla~corbic+flow
#         estfish_bsmt~zoop+flow
#         estfish_bsot~zoop+flow
# '
#2
# modsouth='chla~corbic+flow
#         tzoop~chla+corbic+flow
#         estfish_bsmt~tzoop+flow+corbic+sside
#         estfish_bsot~tzoop+flow+corbic+sside
# '
#3
# modsouth='chla~corbic+flow+temp+turbid
#         tzoop~chla+corbic+flow+temp+turbid
#         amphi~chla+corbic+flow+temp+turbid
#         estfish_bsmt~tzoop+amphi+flow+temp+turbid+corbic+sside
#         #estfish_bsot~tzoop+amphi+flow+temp+turbid+corbic+sside
#         amphi~~tzoop
# '
#4
# modsouth='chla~corbic+flow+temp+turbid
#         hzoop~chla+corbic+flow+temp+turbid
#         pzoop~chla+corbic+flow+temp+turbid+hzoop
#         amphi~chla+corbic+flow+temp+turbid
#         estfish_bsmt~hzoop+pzoop+amphi+flow+temp+turbid+corbic+sside
#         amphi~~hzoop+pzoop
# '
#5
modsouth='chla~corbic+flow+temp+turbid
        hzoop~chla+corbic+flow+temp+turbid
        pzoop~chla+corbic+flow+temp+turbid+hzoop
        fish~hzoop+pzoop+corbic+flow+temp+turbid
        fish=~estfish+estfish_stn+estfish_bsmt
'
modsouth_dtr='chla~corbic+flow+temp+turbid
        hzoop~chla+corbic+flow+temp+turbid
        pzoop~chla+corbic+flow+temp+turbid+hzoop
        fish~hzoop+pzoop+corbic+flow+temp+turbid
        fish=~estfish #+estfish_stn+estfish_bsmt
'
modfitsouth=sem(modsouth, data=filter(fdr,region=="South"))
modfitsouth_dtr=sem(modsouth_dtr, data=filter(fdr_dtr,region=="South"))
summary(modfitsouth, standardized=T, rsq=T)
## lavaan 0.6-10 ended normally after 28 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        30
##                                                       
##                                                   Used       Total
##   Number of observations                            40          47
##                                                                   
## Model Test User Model:
##                                                       
##   Test statistic                                17.243
##   Degrees of freedom                                15
##   P-value (Chi-square)                           0.305
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fish =~                                                               
##     estfish           1.000                               0.912    0.920
##     estfish_stn       0.776    0.124    6.236    0.000    0.708    0.763
##     estfish_bsmt      0.635    0.152    4.170    0.000    0.579    0.587
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   chla ~                                                                
##     corbic            0.284    0.149    1.906    0.057    0.284    0.283
##     flow             -0.124    0.137   -0.906    0.365   -0.124   -0.128
##     temp              0.081    0.140    0.578    0.563    0.081    0.081
##     turbid            0.392    0.158    2.478    0.013    0.392    0.373
##   hzoop ~                                                               
##     chla              0.454    0.112    4.058    0.000    0.454    0.514
##     corbic            0.324    0.110    2.944    0.003    0.324    0.366
##     flow             -0.015    0.098   -0.150    0.881   -0.015   -0.017
##     temp             -0.057    0.100   -0.573    0.566   -0.057   -0.065
##     turbid           -0.028    0.120   -0.236    0.814   -0.028   -0.030
##   pzoop ~                                                               
##     chla              0.510    0.126    4.039    0.000    0.510    0.533
##     corbic           -0.024    0.115   -0.206    0.837   -0.024   -0.025
##     flow             -0.078    0.093   -0.841    0.401   -0.078   -0.084
##     temp              0.056    0.095    0.588    0.556    0.056    0.058
##     turbid            0.180    0.114    1.576    0.115    0.180    0.179
##     hzoop             0.271    0.150    1.802    0.072    0.271    0.250
##   fish ~                                                                
##     hzoop             0.230    0.125    1.833    0.067    0.252    0.222
##     pzoop            -0.095    0.112   -0.849    0.396   -0.104   -0.099
##     corbic            0.174    0.097    1.799    0.072    0.191    0.190
##     flow             -0.080    0.078   -1.021    0.307   -0.088   -0.091
##     temp              0.023    0.080    0.284    0.776    0.025    0.025
##     turbid            0.761    0.102    7.495    0.000    0.835    0.793
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .estfish           0.152    0.069    2.193    0.028    0.152    0.154
##    .estfish_stn       0.360    0.090    3.993    0.000    0.360    0.418
##    .estfish_bsmt      0.640    0.149    4.296    0.000    0.640    0.656
##    .chla              0.723    0.162    4.472    0.000    0.723    0.724
##    .hzoop             0.362    0.081    4.472    0.000    0.362    0.465
##    .pzoop             0.327    0.073    4.472    0.000    0.327    0.358
##    .fish              0.116    0.066    1.777    0.076    0.140    0.140
## 
## R-Square:
##                    Estimate
##     estfish           0.846
##     estfish_stn       0.582
##     estfish_bsmt      0.344
##     chla              0.276
##     hzoop             0.535
##     pzoop             0.642
##     fish              0.860
summary(modfitsouth_dtr, standardized=T, rsq=T)
## lavaan 0.6-10 ended normally after 28 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        25
##                                                       
##                                                   Used       Total
##   Number of observations                            41          47
##                                                                   
## Model Test User Model:
##                                                       
##   Test statistic                                 0.301
##   Degrees of freedom                                 1
##   P-value (Chi-square)                           0.583
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fish =~                                                               
##     estfish           1.000                               0.909    1.000
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   chla ~                                                                
##     corbic            0.120    0.157    0.764    0.445    0.120    0.120
##     flow             -0.033    0.162   -0.202    0.840   -0.033   -0.034
##     temp              0.176    0.158    1.110    0.267    0.176    0.176
##     turbid           -0.048    0.163   -0.291    0.771   -0.048   -0.048
##   hzoop ~                                                               
##     chla              0.443    0.125    3.550    0.000    0.443    0.456
##     corbic            0.279    0.126    2.212    0.027    0.279    0.288
##     flow              0.030    0.130    0.234    0.815    0.030    0.032
##     temp             -0.016    0.128   -0.129    0.898   -0.016   -0.017
##     turbid           -0.148    0.131   -1.137    0.256   -0.148   -0.153
##   pzoop ~                                                               
##     chla              0.543    0.137    3.959    0.000    0.543    0.523
##     corbic           -0.120    0.129   -0.935    0.350   -0.120   -0.116
##     flow             -0.021    0.125   -0.168    0.867   -0.021   -0.021
##     temp              0.162    0.123    1.316    0.188    0.162    0.156
##     turbid           -0.103    0.128   -0.808    0.419   -0.103   -0.099
##     hzoop             0.218    0.150    1.450    0.147    0.218    0.204
##   fish ~                                                                
##     hzoop             0.226    0.157    1.440    0.150    0.248    0.238
##     pzoop            -0.259    0.141   -1.834    0.067   -0.285   -0.292
##     corbic            0.107    0.138    0.778    0.436    0.118    0.117
##     flow             -0.004    0.132   -0.032    0.975   -0.005   -0.005
##     temp              0.120    0.134    0.897    0.370    0.132    0.130
##     turbid            0.348    0.136    2.566    0.010    0.383    0.379
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .estfish           0.000                               0.000    0.000
##    .chla              0.927    0.205    4.528    0.000    0.927    0.950
##    .hzoop             0.591    0.131    4.528    0.000    0.591    0.642
##    .pzoop             0.547    0.121    4.528    0.000    0.547    0.520
##    .fish              0.616    0.136    4.528    0.000    0.747    0.747
## 
## R-Square:
##                    Estimate
##     estfish           1.000
##     chla              0.050
##     hzoop             0.358
##     pzoop             0.480
##     fish              0.253
labelssouth <- createLabels(modfitsouth, cnames)

# residuals(modfitsouth)
# modificationindices(modfitsouth)

Nice plots

Original units

West

North

South

Detrended

West

North

South

Save updated SEM diagrams

plot_grobs <- map(list(figW, figN, figS), ~convert_html_to_grob(.x, 2000))
combined_figure <- ggarrange(plotlist=plot_grobs, labels="auto", 
                             font.label=list(size=9)) %>%
  annotate_figure(top = text_grob("Annual SEM (regional)",
                                  color = "black",
                                  face = "bold",
                                  size = 9))# + 
  #theme(plot.margin = unit(c(1,0,0,0), "lines"))

ggsave('./fig_output/sem_annual_regions.png', combined_figure, width=6, height=6,
       dpi=300)